Understanding fragility in supercooled Lennard-Jones mixtures. 
I. Locally preferred structures 
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We reveal the existence of systematic variations of isobaric fragility in different supercooled 
Lennard-Jones binary mixtures by performing molecular dynamics simulations. The connection 
between fragility and local structures in the bulk is analyzed by means of a Voronoi construction. 
We find that clusters of particles belonging to locally preferred structures form slow, long-lived do- 
mains, whose spatial extension increases by decreasing temperature. As a general rule, a more rapid 
growth, upon supercooling, of such domains is associated to a more pronounced super- Arrhenius 
behavior, hence to a larger fragility. 

PACS numbers: 61.43.Fs, 61.20.Lc, 64.70.Pf, 61.20.Ja 



I. INTRODUCTION 

What is the origin of super- Arrhenius behavior of 
transport coefficients and relaxation times in supercooled 
liquids? This is one of the long-standing, open ques- 
tions regarding the physics of the glass-transition. An- 
gelic introduced the notion of fragility to quantify the 
degree of super-Arrhenius behavior: generally speaking, 
the steeper the increase of relaxation times by decreasing 
temperature, the more fragile the glass-former. While re- 
laxation times increase by several orders of magnitude on 
approaching the glass-transition, only mild variations in 
the average liquid structure are observed. Nevertheless, a 
subtle, deep link between microstructural order and dy- 
namics is believed to existi^'^'^ In particular, the emer- 
gence, upon supercooling, of slow domains, characterized 
by well-defined locally preferred structures has been re- 
cently identified as a possible origin of super-Arrhenius 
behavior, both in numerical simulations^'^ and theoret- 
ical approaches. ^''^'^^ Despite some efforts in this direc- 
tion, the microscopic foundations of such a connection 
have been explored only to a limited extent. In fact, 
the relation between fragility and local order has been 
explicitly investigated only for monoatomic systems^^ii 
and operational schemes for a direct determination of lo- 
cally preferred structures have not been tested yet for 
model supercooled liquids.-^'-^ 

Understanding how fragility changes in different sys- 
tems appears to be an even harder task. Experi- 
mentally, some interesting correlations have been pro- 
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but a sharp interpretation of these 



results is hindered by the complexity and the varying 
nature of intermolecular interactions. Numerical sim- 
ulations, allowing fine-tuning of interaction parameters 
for a wide choice of potentials, are optimal tools for 
studying how the details of the interaction affect the be- 
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havior of supercooled liquids. Following this approach, 
some authors have recently investigated the connection 
between the features of the interaction potential and 
fragility.— i2Si2iiS^ For instance, it has been shown that 
changing the power-law dependence of repulsion in a soft- 
sphere mixture leaves fragility invariant while altering 
both repulsive and attractive parts of the Lennard-Jones 
potential can change the fragility of a supercooled mix- 
ture,^'' The role of intermolecular interactions should not, 
however, be overemphasized. Statistical mechanics theo- 
ries of the glass-transition often use correlation functions 
or other coarse-grained information on the liquid proper- 
ties as input, rather than the bare interaction potential. 
For instance, the most celebrated mode coupling theorjs^^ 
emphasizes the role of pair correlations, as described by 
the static structure factors. Providing correlations be- 
tween fragility and properties such as local order, which 
may be used to characterize mesoscopic domains in su- 
percooled liquids, would also be more helpful for under- 
standing experimental trends. 

In this paper, we investigate the connection between 
fragility and local order by performing Molecular Dynam- 
ics simulations of different supercooled Lennard-Jones 
mixtures. We consider a set of equimolar, additive mix- 
tures with varying size ratio, together with some proto- 
typical mixtures (Sec.[TT|- Different from most numerical 
simulations on Lennard-Jones mixtures, we cool these 
systems at constant pressure (Sec. Will . The isobaric 
fragilities obtained from relaxation times and diffusion 
coefficient show a systematic variation in additive mix- 
tures upon varying size ratio (Sec. lIV[) . Different fragility 
indexes are also found for two models2ii2^ that have of- 
ten been employed in computational studies on the glass- 
transition. We rationalize these findings by analyzing 
the role of locally preferred structures, as identified by 
means of a Voronoi construction (Sec. |V|. We show that 
super-Arrhenius behavior of dynamical properties can be 
ascribed to a rapid growth, upon supercooling, of slow 
domains possessing distinct microstructural features. As 
a keynote, we find that the fragility of binary mixtures 
is correlated to the thermal rate of growth of such do- 
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mains: the more fragile the mixture, the more rapid the 
increase of the fraction of particles forming locally pre- 
ferred structures. 



II. MODEL BINARY MIXTURES 

We have performed extensive Molecular Dynamics sim- 
ulations for 13 binary Lennard-Jones mixtures. All mod- 
els are composed of = 500 classical particles enclosed 
in a cubic box with periodic boundary conditions. Par- 
ticles interact via the Lcnnard- Jones potential 



(1) 



where a,P = 1,2 are indexes of species. In our conven- 
tion, particles of species 2 have a smaller diameter than 
those of species 1 {a22 < fn), and we fix an = 1.0 for 
all systems. Reduced units will be used in the following, 
assuming an, en and ^/rnjof^Jen as units of distance, 
energy and time respectively. We have employed the cut- 
off scheme of Stoddard and Fordj^- which ensures con- 
tinuity up to the first derivative of Ua/3{r) at the cutoff 
radius Tc — 2.5. 

The interaction parameters of the mixtures considered 
in this work are shown in Table IH These models are 
characterized by a varying degree of non-additivity in 
the composition rule for the cross interaction and by dif- 
ferent values of size ratio A — a^ilow. Deviations from 
Lorentz-Berthelot composition rules can be quantified us- 
ing 
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We will now proceed to a brief presentation of the rele- 
vant features of these force fields. 



BMLJ — Binary Mixture of Lennard-Jones particles 

This is the classic mixture of Kob and Andersen,^** 
which has been used extensively as a model supercooled 
liquid. It is characterized by a significant asymmetry 
in the interaction parameters, both in the interaction 
diameters and in the energy scales of the two species. 
Furthermore, cross interactions are strongly non-additive 
(77 = 0.85, 77 — 2.1). The number concentration of large 
particles is fixed at x\ = 0.8, as in the original work. 



TABLE L Parameters of Lennard-Jones potentials for binary 
mixtures. Also shown are the masses m\ and 7712 of the two 
species and the concentration x\ of particles of species 1. In 
the case of additive mixtures AMLJ-A, the following values of 
size ratio A have been used: 0.60, 0.64, 0.70, 0.73, 0.76, 0.82, 
0.88, 0.92, 0.96, 1.00. 
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a single energy scale is present. The masses of the two 
chemical species are equal. The number concentration 
X\ = 0.67 allows deep supercooling of the mixture. 



WAHN — Additive mixture of Wahnstrom 

Introduced by Wahnstrom, this equimolar mixture 
is actually the Lennard-Jones version of the super- 
cooled soft-sphere model used in the early simulations 
of Bcrnu et al.?^ It has been employed several times in 
the literature as a model glass-former. The interaction 
parameters are additive and characterized by a moderate 
size asymmetry (A — 0.837). Note that the mass ratio is 

TflilTflY = 0.5. 



AMLJ-A — Additive Mixture of Lennard-Jones 
particles 

This is a set of equimolar additive mixtures. The 
masses of the two species are equal m\ = TO2 = 1.0 and 
the size ratio A is allowed to vary, keeping cn fixed at 
1.0. Recently, Lennard-Jones clusters of this type has 
been investigatedi^ Ten different values of A have been 
used in the range 0.60 < A < 1.0. Note that AMLJ-0.837 
would be the same as the WAHN mixture, if not for the 
the different mass ratio. 



NiasYe?— Lennard-Jones model for NissYer alloy 



III. QUENCHING PROTOCOLS AND 
SIMULATION DETAILS 



The parametrization of this mixture has been intro- 
duced by Delia Valle et al^ to provide a realistic de- 
scription of the structural features of binary amorphous 
alloys of Ni and Y atoms. The cross-interaction diameter 
is non-additive (77 = 0.91), as in the case of BMLJ, but 



Numerical simulations of model supercooled liquids are 
usually performed at constant density. Recently, some 
attention has been drawn to density effects^i^ and to 
the role of pressure i^i^ Nonetheless, numerical stud- 
ies of supercooled Lennard-Jones mixtures along isobaric 
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FIG. 1: Bathia-Thornton structure factors in the deeply supercooled regime. Number-number (solid lines), number- 
concentration (dashed lines) and concentration-concentration (dotted lines) structure factors are shown. The concentration- 
concentration structure factor has been normalized to one by plotting Sec (k) /{xiX2). Data are shown at the lowest equilibrated 
temperature for each given system. Note that all structure factors are finite in the limit A: ^ and that the first sharp peak 
of S'jvjv(fc) around feo ~ 8 is roughly system independent. 
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FIG. 2: Temperature dependence of density p{T) along iso- 
baric quenches at P = 10 (left axis, circles) and of pressure 
P{T) in isochoric quenches (right axis, squares). Data are 
shown for BMLJ (filled symbols) and WAHN (open symbols) . 
Isochoric quenches were performed at p = 1.2 for BMLJ and 
p = 1.3 for WAHN. 



quenches have seldom been performed. In our case, fol- 
lowing cooling paths at a common, constant pressure pro- 
vides a simple mean to perform a homogeneous compar- 



ison of different mixtures. Moreover, this approach is 
closer to the one usually employed in experimental condi- 
tions, where most comparisons of supercooled liquids are 
performed at constant (atmospheric) pressure. Whereas 
some minor differences between isochoric and isobaric 
quenches can be observed, and will be highlighted in the 
following, we remark that the main conclusions of this 
work are independent of the choice of the quenching pro- 
tocol. 

Isobaric quenches were performed by coupling the sys- 
tem to both Berendsen thermostat and barostat during 
the equilibration phase. In most cases, we employed the 
standard Velocity- Ver let algorithm for production runs. 
In order to achieve better control on temperature in the 
deeply supercooled regime, we also performed a few pro- 
duction runs using the Nose-Poincare thermostat i^i^^ 
For the system size and temperatures in consideration, 
simulations in different ensembles provide consistent re- 
sults, as far as average dynamical and thermodynamical 
properties are concerned. Isochoric quenches were per- 
formed in a similar way, using the Berendsen thermo- 
stat only during equilibration. The timestep St of both 
Velocity- Verlet and Nose-Poincare integrators was varied 
from 0.002 at high temperature to 0.006 in the super- 
cooled regime. The time constant for the Berendsen ther- 
mostat'^^ was tx = St/ 0.1, while the coupling constant 
the for Berendsen barostat^ was lO'^ in reduced units. 
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The inertia parameter of the Nose-Poincare thermostat 
was set to Q = 5. Thanks to the symplectic natm'e of 
the integrators and to the smooth cutoff employed, no 
significant drift in the relevant conserved Hamiltonian is 
observed in either kind of simulation used for production 
runs, even for very long runs (up to 2.5 x 10^ steps at the 
lowest equilibrated temperatures). 

Beside checking the stability of mean potential energy 
and pressure during the production phase, we adjusted 
the total simulation times in order to achieve similar val- 
ues of total root mean square displacement for all state 
points and for all mixtures. Typical values of total root 
mean square displacement for large particles are always 
larger than Aau. To check the reliability of our results, 
we also tried both faster and slower cooling rates, and 
performed some thermal histories by reheating deeply su- 
percooled samples. In the case of AMLJ-A mixtures, we 
found that they could be safely supercooled for values of 
size ratio A between 0.60 and 0.84. In the case A = 0.88, 
in fact, some early signs of crystallization were found 
in our sample. A part from this case, no sign of phase 
separation or crystallization was detected looking at the 
time evolution of thermodynamic and structural proper- 
ties. The Bathia-ThorntonSI number-number, number- 
concentration and concentration-concentration structure 
factors for a selection of mixtures are shown in Fig. [T] at 
the lowest equilibrated temperatures. 

In the following, we will mostly consider isobaric 
quenches performed at a reduced pressure P ~ 10. 
Reference data for BMLJ are available along this iso- 
barjSS so that we could check the reliability of both 
dynamical and thermodynamical properties obtained in 
our simulations. In the case of BMLJ, we performed 
a few isobaric quenches over a wider range of pressure 
(P = 5, 10, 20, 50) to investigate the pressure dependence 
of isobaric fragility. Isochoric quenches have been per- 
formed for BMLJ and WAHN mixtures, fixing the den- 
sity at p = 1.2 and p = 1.3 respectively. These values of 
p are equal to the ones used in the original papersi2ii2^ 
In Fig. [5] we show the temperature dependence of density 
(at constant pressure) and pressure (at constant density) 
for BMLJ and WAHN. Note that the range of pressure 
investigated for BMLJ (5 < P < 50) is consistent with 
the variation of pressure for this system along the iso- 
chore p — 1.2. 



IV. FRAGILITY 

The fragility of a glass-former quantifies how rapid the 
change of dynamical properties is upon supercooling. We 
will focus on dynamical quantities that can be computed 
with good statistical accuracy in numerical simulations, 
namely relaxation times for the decay of density fluctua- 
tions and diffusion coefficients. 

For the definition of relaxation times we will consider 



the self part of the intermediate scattering function 
Kik, i) = ]^ E< • [^^(^ + - ) (3) 

^ i— 1 

where a = 1, 2 is an index of species and ( ) indicates 
an average over time origins Iq. Relaxation times for 
species a are defined by the condition F"{k*,Ta) = 1/e, 
where k* corresponds to the position of the first peak in 
the number- number structure factor (see Fig. [1]). The 
value of k* is roughly system- and temperature indepen- 
dent for the mixtures in consideration, and close to fc « 8. 
In the following, we will focus on the temperature depen- 
dence of T = Ti, but similar trends are observed when 
considering the small particles. 

The difficulty of providing an unbiased, global descrip- 
tion of the temperature dependence of transport coeffi- 
cients and relaxation times by fitting the experimental 
data has been particularly stressed by Kivelson et al.t^ 
Care must be taken when the definition of fragility itself 
relies on a specific functional form, or when the latter 
is used for extrapolations outside the accessible range of 
temperature. We thus seek functional forms that are re- 
liable over a large range of temperature and require the 
range for fitting to be well-specified and physically mo- 
tivated. For describing the temperature dependence of 
relaxation times, we start with the well-known Vogel- 
Fulcher-Tammann (VFT) law and write it in the form^ 



r(T) 



> exp 



1 



K{T/To-l) 



(4) 



The material-dependent parameter K quantifies the 
fragility of the glass-former under consideration. The 
larger is K, the steeper is the increase of t{T) upon su- 
percooling. Equation ^ provides a fairly good descrip- 
tion of relaxation times in the deeply supercooled regime, 
but is inaccurate at high temperature.'^® In the normal 
liquid regime, in fact, relaxation times have a mild tem- 
perature dependence, which is well described by the Ar- 
rhenius law. The existence of a crossover between these 
two regimes around some temperature Tonset, accompa- 
nied by several qualitative changes in the properties of 
the liquid, is well established in the literature^®. It seems 
thus sensible to use the following global functional form 



r(T) = 



Too exp 
r^exp 



Eoc/T] 

1 



K{T/To-l) 



where 



, exp 



Eoo/T* 



1 



K{T*/To - 1) 



(5) 



(6) 



as a generalized VFT law. This functional form is con- 
tinuous at the crossover temperature T* and provides 
a fragility index K with analogous physical meaning to 
that in Eq. (jH). To fit our simulation data to Eq. (g]). 
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TABLE 11: Fitted parameters for relaxation times t{T) of large particles according to the generalized Vogel-Fulcher-Tammann 
given by Eq. ([Sj, and for total diffusion coefficient D(T) according to Eq. The reference temperature T,- and the onset 
temperature Tonaet are described in the text. 



Relaxation times Diffusion coefficient 





P 


Ton set 


Tr 


Too 


Eoc 


To 


K 


To 


K 


BMLJ 


5.0 


0.75 


0.464 


0.093f(3) 


f.99(f) 


0.392(3) 


0.43(2) 


0.36f(6) 


0.37(2) 


BMLJ 


fO.O 


0.95 


0.574 


0.08f5(5) 


2.6f(f) 


0.479(2) 


0.40(f) 


0.457(3) 


0.4f(l) 


BMLJ 


20.0 


f.20 


0.765 


0.067(f) 


3.7f(9) 


0.63(f) 


0.40(5) 


0.60(f) 


0.42(3) 


BMLJ 


50.0 


f.80 


f.248 


0.048f(9) 


6.60(7) 


f.03(f) 


0.38(3) 


f.Of(f) 


0.46(2) 


WAHN 


fO.O 


0.90 


0.623 


0.0825(4) 


2.33(f) 


0.573(4) 


0.94(6) 


0.523(6) 


0.70(4) 


WAHN 


20.0 


f.20 


0.825 


0.0670(6) 


3.38(3) 


0.752(5) 


0.84(5) 


0.697(9) 


0.68(4) 


NiasYeT 


fO.O 


0.90 


0.489 


0.0777(7) 


2.53(2) 


0.39f(3) 


0.3f(f) 


0.379(4) 


0.37(1) 


AMLJ-0.60 


fO.O 


0.85 


0.45f 


0.076(f) 


2.43(3) 


0.34f(5) 


0.24(f) 


0.319(4) 


0.29(1) 


AMLJ-0.64 


fO.O 


0.90 


0.474 


0.0769f(3) 


2.444(f) 


0.38f(4) 


0.32(f) 


0.356(4) 


0.36(1) 


AMLJ-0.70 


fO.O 


0.90 


0.5f4 


0.08ff(f) 


2.359(4) 


0.440(5) 


0.48(3) 


0.417(7) 


0.51(4) 


AMLJ-0.73 


fO.O 


0.90 


0.560 


0.0785(7) 


2.48(2) 


0.502(4) 


0.7f(4) 


0.466(6) 


0.62(4) 


AMLJ-0.76 


fO.O 


0.90 


0.60f 


0.0790(9) 


2.49(3) 


0.554(4) 


f.0f(9) 


0.519(7) 


0.84(7) 


AMLJ-0.82 


fO.O 


0.95 


0.636 


0.0803(5) 


2.53(f) 


0.59f(6) 


1-1(1) 


0.53(f) 


0.76(7) 




P 


Ton set 


Tr 






To 


K 


To 


K 


BMLJ 


1.2 


f.OO 


0.422 


0.ff0(2) 


2.69(2) 


0.33f(3) 


0.3f(f) 


0.328(2) 


0.426(9) 


WAHN 


f.3 


1.05 


0.522 


0.097(f) 


2.73(2) 


0.447(5) 


0.50(4) 


0.426(2) 


0.64(1) 



we proceed similarly to Kivelson et al..^^ First we fit the 
relaxation times to the Arrhenius law Too exp(£'(x;/r) in 
the range T > Tonset and then we use t^o and E^o as 
fixed parameters in a global fit to Eq. (O. In this way, 
a good overall fit is obtained for relaxation times. Note 
that T* is considered as a fitting parameter, but the main 
conclusions of this section will not be altered by fixing 

— onset • 

The Angell plots of relaxation times in Fig.[3]and Fig. [3] 
constitute the starting point of our discussion. For their 
construction we have used a reference temperature T^, 
akin to the glass transition temperature Tg, at which 
t(Tj.) = 4 X 10'^. Such a value of T{Tr) is close to the one 
used by Bordat et al^ in a study of modified Lennard- 
Jones mixtures. For each mixture, the value of the ref- 
erence temperature T^, which we have extrapolated us- 
ing Eq. (O, is only slightly below the lowest equilibrated 
temperature. 

In Fig.[3]we consider the set of additive, equimolar mix- 
tures AMLJ-A along isobaric quenches at P = 10. The 
size ratio A is varied in the range 0.60 < A < 0.82. A 
strong, systematic variation is apparent upon varying A: 
the mixture becomes more fragile as A increases, i.e. as 
the size asymmetry between the two species is reduced. 
Recently, a similar influence of size ratio on fragility has 
been observed in modified BMLJ mixtures. The trend 
of variation of fragility is confirmed by our fitting pro- 
cedure, whose outcome is summarized in Table |lTl The 
dependence of the isobaric fragility index K on A, shown 
in the inset of Fig. [31 also suggests the existence of a sat- 
uration of fragility around A = 0.80. This feature will be 
further discussed in Sec. |Vl in connection with icosahe- 
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FIG. 3: Angell plot of relaxation times of large particles r for 
a selection of AMLJ-A mixtures. Results are shown for A = 
0.60, 0.70, 0.73, 0.82 along the isobar P = f 0. The reference 
temperature Tr is described in the text. The inset shows 
the isobaric fragility index K obtained from generalized VFT 
equation (see Eq. ([5}) against size ratio A. 



dral ordering. 

Figure [4] shows the temperature dependence of relax- 
ation times for the two well-studied glass-formers BMLJ 
and WAHN, both cooled isobarically at P = 10. These 
mixtures have been used extensively for numerical inves- 
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FIG. 4: Angell plot of relaxation times of large particles r 
for BMLJ (black circles) and WAHN mixture (white circles) 
along isobaric quenches at P = 10. The inset shows results 
at P = 5, 10, 20, 50 for BMLJ. 



FIG. 5: Angell plot of total diffusion coefficient D for a se- 
lection of AMLJ-A mixtures. Results are shown for A = 
0.60,0.70,0.73,0.82 along the isobar P = 10. The reference 
temperature To is obtained from fit to Eq. ((Tjl. The inset 
shows the isobaric fragility index K obtained from Eq. ([7]) 
against size ratio A. 



tigations of the glass-transition, but a direct comparison 
has never appeared in the literature. We find that WAHN 
is appreciably more fragile than BMLJ, independent of 
quenching protocols and system sizeSS,. The enhanced 
fragility of WAHN is not surprising since this mixture is, 
a part from a different mass ratio, an AMLJ-A mixture 
with A — 0.837. We found that such difference in mass 
ratio is irrelevant to the dynamical properties in consider- 
ation. Thus, WAHN can be considered as the end-point 
of a series of supercooled mixtures with increasingly large 
fragility. 

Does isobaric fragility itself depend on pressure? 
This question has received much attention in the last 
years, in particular within the experimental commu- 
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A tentative answer can be given for 



BMLJ, for which we performed isobaric quenches in the 
range 5 < P < 50. By looking at the inset of Fig. 21 
we can see that relaxation times obtained along different 
isobars collapse on a master curve by scaling T with the 
corresponding Tr. Numerical values of isobaric fragility 
K , obtained from Eq. ([5]) along different isobars, are also 
very close to each other (see Table lll| . Thus, our results 
indicate that the pressure dependence of isobaric fragility 
of Lennard- Jones mixtures might be mild or even negli- 
gible for moderate variations of P. For a given system, 
we also find that the isochoric fragility is slightly smaller 
than the corresponding isobaric fragility, in agreement 
with experimental observations4^ 

As another dynamical indicator we consider the the 
total diffusion coefficient D = xiDi +X2D2, given by the 
usual Einstein relation. The temperature dependence of 
D can be described quite satisfactorily over the entire 
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FIG. 6: Angell plot of total diffusion coefficient D for BMLJ 
(black circles) and WAHN mixture (white circles) along iso- 
baric quenches at P = 10. The inset shows results at 
P = 5, 10, 20, 50 for BMLJ. 



temperature range by a modified VET law^ 

1 1 



D{T) = DoTexp 



K{T/To - 1) 



(7) 



Again, the parameter K in Eq. ([7]) provides a measure 
of fragility. Angell plots for the diffusion coefficient are 
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shown in Figs. [5]and[6l and results of fits to Eq. ([7]) are 
collected in Table HIl Re garding the variation of fragility 
in our models, results obtained for the diffusion coef- 
ficient confirm the analysis based on relaxation times. 
Note that, by considering the temperature dependence 
of both r = Ti and D, we account for a possible depen- 
dence of fragility on chemical species and fc- vector. Thus, 
the following trends appear to be rather robust: (i) The 
increase of fragility with size ratio A in AMLJ-A mix- 
tures, (ii) The less fragile behavior of BMLJ compared 
to WAHN. (iii) The independence of isobaric fragility on 
pressure in BMLJ. 



V. LOCALLY PREFERRED STRUCTURES 

The growth of long-lived, slow domains, characterized 
by well-defined structural features, has been discussed re- 
cently in connection to super-Arrhenius behavior of dy- 
namical properties of supercooled liquids i^i^i^ii^ In this 
section, we will follow a similar approach by analyzing 
the properties of local order in the Lennard- Jones mix- 
tures introduced in Sec. |TT1 We will highlight the exis- 
tence of different, well-defined local geometries in these 
systems and show how their properties are related to the 
variations of fragility. 

For this kind of analysis we have employed a Voronoi 
construction.""^ Each particle in the system is the cen- 
ter of a Voronoi polyhedron, which is constructed by 
intersection of planes orthogonal to all segments con- 
necting the central particle to the other ones. Planes 
are drawn at a fraction faf3 of these segments, where 
a is the species of the central particle and (3 is the 
species of the other particle. We have used the recipe 
/q/3 = o'aa/ if^aa + o'/i/s),^ but the main conclusions of 
this section will not be altered when considering the more 
intuitive choice, fafs = 1/2. The sequence (n3,n4,...), 
where Uk is the number of faces of the polyhedron hav- 
ing k vertices, provides a detailed description of the local 
geometry around a given particle. All zero values behind 
the maximum number of vertices of a polyhedron are ig- 
nored. We have applied this kind of analysis to both in- 
stantaneous configurations, sampled along Molecular Dy- 
namics trajectories, and to local minima of the potential 
energy, obtained by conjugate-gradients minimizations. 
Between 200 and 2000 independent configurations have 
been analyzed for each state point. We found that it is 
easier to characterize local order in the Lennard-Jones 
mixtures considered in this work by taking "the point of 
view" of small particles. Well-defined geometries appear, 
in fact, most frequently around small particles, while no 
recognizable local order is apparent around large parti- 
cles. In the following, wc will thus concentrate our at- 
tention on the properties of Voronoi polyhedra centered 
around small particles. 

Such a Voronoi construction will provide an indica- 
tion of what are the locally preferred structures^iii of 
our model supercooled liquids. In the rest, we will use 



TABLE III: Most frequent Voronoi polyhedra around small 
particles. The percentage is computed with respect to the 
number of small particles in the system. Also shown is the 
average number of neighbors of species 1 (ni) and 2 {112) ■ Re- 
sults refer to local minima along the isobar P = 10 and are 
shown for T — 2.0 and slightly above the reference tempera- 
ture Tr, i.e. for the lowest equilibrated temperature. 
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an effective, purely geometric definition of locally pre- 
ferred structures, as the ones corresponding to the most 
frequent polyhedra found in our Voronoi construction. 
Determining unambiguously the origin of the preference 
for a given local structure in the bulk will require a signif- 
icant additional effort, since, in general, such preference 
will depend in a non-trivial way on the environment sur- 
rounding a given local structure and may be triggered 
by factors other than energetic stability. For instance, 
packing effects can play an important role in stabilizing a 
local structure. Moreover, in the case of binary mixtures, 
compositional freedom further increases the complexity 
of such analysis. Taking into account the effect of the liq- 
uid environment around a local structure at a mean-field 
level, and identifying the appropriate "local free energy" 
to be minimized, remains an open issue of current re- 
searchJii^ Here, we tackle these issues by using a purely 
geometric definition of locally preferred structures, which 
has the further advantage of being available even when 
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FIG. 7: Temperature dependence of the fraction of small par- 
ticles at the center of (0,0,12)-polyhedra in AMLJ-A for se- 
lected values of A. Results are shown for instantaneous config- 
urations (main plot) and local minima (inset) along isobaric 
quenches at P = 10. 



FIG. 8: Variation of icosahedral ordering with size ratio in 
additive mixtures AMLJ-A. The fraction of small particles at 
the center of (0,0,12)-polyhedra in local minima is shown as a 
function of size ratio A, at T = 2.0 (black circles) and at the 
lowest equilibrated temperatures (open circles). 



energetic criteria will fail, e.g. for hard-spheres. 

Let us first focus on equimolar, additive mixtures 
AMLJ-A. One of the most relevant structural features 
of these mixtures is the existence of a varying degree 
of icosahedral ordering. Icosahedral coordination has a 
sharp signature in the Voronoi construction, being associ- 
ated to (0,0,12)-polyhedra, i.e. 12 pentagonal faces. The 
temperature dependence of the fraction of icosahedra for 
different AMLJ-A mixtures, shown in Fig. [3 displays a 
striking correlation with the variation of fragility. In 
fact, the increase of icosahedral ordering upon supercool- 
ing is more rapid and more pronounced as A increases, 
i.e. as the fragility of the mixture increases^iS.. To our 
knowledge, this is the first time that such a relationship 
is established in supercooled binary mixtures. Previous 
numerical studies have focused, in fact, on the connec- 
tion between icosahedral ordering and super- Arrhenius 
behavior in monoatomic liquidsi^ii The theoretical inter- 
pretation of our results is by no means trivial. On the 
one hand, the variation of fragility with icosahedral or- 
dering may be understood within the frustration-limited 
domains theory^ in terms of a more rapid stabilization, 
upon supercooling, of locally preferred structures — in 
the present case, icosahedra — . We will further discuss 
this point below. On the other hand, the trend we find 
in our simulations and shown in Fig. [7] appears to be 
at variance with the results of a recent phenomenologi- 
cal modelj^ii^ Further clarifications on the relevance of 
this theoretical approach to our simulated systems are 
required. 

To have a better feeling of how icosahedral ordering 
is triggered by size ratio, we show, in Fig. [8l the frac- 
tion of (0,0,12)-polyhedra in local minima as a function 
of A. Results are shown along the isotherm T — 2.0 
and for T k Tr, i.e. at the lowest temperatures that 



could be accessed in equilibrium condition. Depending 
on temperature, a different range of A is considered. In 
the deeply supercooled regime (T « Tj.) only mixtures 
with 0.60 < A < 0.84 could be equilibrated (see Sec. lnT| . 
For variation of A in this range, icosahedral ordering in- 
creases with size ratio, in a way which strongly resem- 
bles the increase of fragility index K with A, reported in 
Sec. |IV1 At high temperature (T = 2.0) the full range 
0.60 < A < 1.00 can be accessed and our data reveal the 
existence of a maximum of icosahedral ordering around 
A w 0.82. This feature might provide a simple explana- 
tion to the existence of a saturation of fragility around 
A — 0.80 reported in Sec. |TVl Interestingly, the results 
obtained from local minima at high temperature show 
that the variation of icosahedral ordering with size ratio, 
which is apparent in the deeply supercooled regime, is 
already encoded in the liquid inherent structure4i 

Such a pattern of variation of icosahedral ordering 
with size ratio is strikingly similar to that observed in 
models of bidisperse Cu glasses^ and in the realistic 
models of metallic glasses developed by Hausleitner and 
Hafner42.i^ This suggests that the increase of icosahe- 
dral ordering with size ratio, and its consequent correla- 
tion with fragility, might be a general feature of binary 
systems with additive, or nearly additive, spherical in- 
teractions^. Furthermore, the onset of crystallization 
for A > 0.88 could be simply be explained, in this kind 
of systems, by the decrease of icosahedral coordination 
and by a larger occurrence of (0,3,6,4)-polyhedra and 
(0,4,4,6)-polyhedra, which are typical of FCC crystals. 51 
We found, in fact, that the fraction of FCC-related poly- 
hedra in the normal liquid regime increases as A — > 1. A 
similar behavior has been explicitly demonstrated for the 
bidisperse Cu model by means of a Honeycutt- Andersen 
construction. 

We can now generalize the connection between fragility 
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(a) (0,2,8)-polyhedron (b) (0,3,6)-polyhedroii 




(c) (0,0,12)-polyhedroii 



FIG. 9: Examples of locally preferred structures found in lo- 
cal minima of supercooled Lennard- Jones mixtures. Small 
and large particles are shown as dark and pale spheres re- 
spectively, (a) (0,2,8)-polyhedron (twisted bicapped square 
prism) in BMLJ. This is the most frequent chemical coordina- 
tion, incidentally one or two small particles can form the cap. 

(b) (0,3,6)-polyhedron (capped trigonal prism) in NiasYer . In 
this case, one of the caps is often formed by a small particle. 

(c) (0,0,12)-polyhedron (icosahedron) in WAHN. On average, 
the coordination around the central particle is equimolar. 



and locally preferred structures by analyzing systems 
possessing favored geometries different from icosahedra. 
Such opportunity is provided by non-additive mixtures, 
such as BMLJ and NisaYgy . In the case of NisaYgv , in 
fact, icosahedral ordering has been shown to be strongly 
frustratedi^i Non-additivity of the interaction potential 
favors the formation of trigonal prismatic structures f21 
similar to those found in the crystalline phases of NiY 
alloys. This is confirmed by our analysis, which shows 
(see Table IIIip that the most frequent Voronoi polyhe- 
dron around small particles in NissYgTis the (0,3,6)- 
polyhedron, corresponding to capped trigonal prismatic 
structures. In the case of BMLJ, we find that the (0,2,8)- 
polyhedron has the largest occurrence around small par- 
ticles both at high and low temperature (see Table HIT)) . 
This polyhedron corresponds to twisted bicapped square 
prisms, mostly formed by neighboring large particles. 
The preference for twisted prismatic structures in BMLJ 
has also been highlighted in studies on the coordina- 
tion polyhedra in the liquid, and on the stability of 
isolated clusters. As a working hypothesis, we iden- 
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FIG. 10: Bond-angle distributions around small particles for 
WAHN (upper plot), BMLJ (middle plot) and NiasYe? (lower 
plot). The bond-angle distribution fi2i{0) is shown as solid 
line. Also shown is the bond-angle distribution fi2i{0) re- 
stricted to small particles which are at the center of the locally 
preferred structure of the system, as given by Fig. [9] Data 
refer to the lowest equilibrated temperature of each given sys- 
tem. The sharp peaks in the fi2i{0) distributions filtered for 
locally preferred structures reflect the ideal angles of the cor- 
responding geometry. 



tify the geometries associated to (0,3,6)-polyhedra and 
(0,2,8)-polyhedra as the locally preferred structures of 
NiasYgr and BMLJ, respectively. These two non-additive 
mixtures can be effectively contrasted to WAHN, which 
displays a strong icosahedral ordering upon supercooling. 

In Fig. [9] we show three highly symmetric configura- 
tions corresponding to locally preferred structures, found 
in local minima of WAHN, BMLJ and NisaYey . Notice 
that the structures shown in the figure are among the 
most symmetric in their own class of Voronoi polyhedra. 
Support to our definition of locally preferred structures 
is given by the analysis of angular distributions. To this 
aim, we compute the bond-angle distribution functions 
faf3-f{0) around particles of species a, where (3 and 7 are 
the species of the neighboring particles. Particles shar- 
ing a face in the Voronoi construction are considered as 
neighbors. In Fig. [TO] we focus on the bond-angle distri- 
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FIG. 11: Temperature dependence of the fraction of small 
particles at the center of selected Voronoi polyhedra in in- 
stantaneous configurations (main plot) and local minima (in- 
set). The fraction of (0,2,8)-polyhedra in BMLJ (white sym- 
bols) and (0,0,12)-polyhedra in WAHN mixture (black sym- 
bols) are shown along isobaric quenches at P = 10. Data for 
NissYer are close to those for BMLJ, but are not shown for 
clarity. 



bution fi2i{0) for central small particles and neighboring 
large particles in WAHN, BMLJ and NiasYgr . We com- 
pare the total fi2i{d) to the one restricted to small parti- 
cles being the center of the typical Voronoi polyhedron of 
the mixture. A strong resemblance is observed between 
the average environment around small particles in the 
bulk and the local geometries of Voronoi polyhedra cor- 
responding to our putative locally preferred structures. 

As it can be seen from Fig. [Til the thermal rate of 
growth of the fraction of particles forming locally pre- 
ferred structures is again correlated to the fragility of 
the model. In fact, the fraction of icosahedra in WAHN 
increases rapidly by decreasing temperature, whereas the 
growth of prismatic structures, typical of BMLJ and 
NiasYgT , is rather mild. The most frequent Voronoi 
polyhedra of all these mixtures are not homogeneously 
spread in the system. They tend to form growing do- 
mains as the temperature decreases. This feature is ex- 
emplified by the two snapshots in Fig. ll2[ where we show 
the typical extension of domains formed by locally pre- 
ferred structures in local minima, for deeply supercooled 
BMLJ and WAHN mixtures. Similar extended domains 
formed by interlocking icosahedra have been found in the 
supercooled regime of Dzugutov liquids.— ^ 

To assess the statistical relevance of the presence of 
such domains, we analyze the distribution P{N) of clus- 
ters composed by N neighboring particles forming locally 
preferred structures. Our identification of domains is as 
follows. For each given configuration, we partition the 
particles into three classes: (i) c-particles, which are the 
center of a locally preferred structure; (ii) n-particles, 
which are neighbors to some other c-particle, but are not 
themselves centers of a locally preferred structure; (iii) 




(a) BMLJ (P = 10, T = 0.60) 




(b) WAHN (P = 10, r = 0.645) 

FIG. 12: Domains formed by locally preferred structures in 
local minima at the lowest equilibrated temperature at P = 
10 (WAHN: T = 0.645. BMLJ: T = 0.60). Particles forming 
(a) (0,2,8)-polyhedra in BMLJ and (b) (0,0,12)-polyhedra in 
WAHN are shown as spheres of the same radius, irrespectively 
of chemical species. 



o-particles, which are neither c-particles nor n-particles, 
i.e. they are outside the domains formed by locally pre- 
ferred structures. Particles sharing a face in the Voronoi 
construction are then considered as neighbors. Using the 
partitioning scheme above, we identify domains as clus- 
ters composed by neighboring c- and n- particles. The 
distribution P{N) is shown in Fig. [13] for instantaneous 
configurations of BMLJ and WAHN at different state 
points. By decreasing temperature, a clear tendency of 
forming larger clusters is observed in both systems. In 
the WAHN mixture we find that, around T^., there is 
almost always a large cluster formed by icosahedra per- 
colating in the simulation box, beside some smaller ones. 
This feature is reflected in the bimodal distribution of 
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FIG. 13: Distribution P{N) of the size A'' of domains formed by locally preferred structures in BMLJ (left plots) and WAHN 
(right plots) at different T. Results refer to isobaric quenches at P = 10. 



P{N) for WAHN at low temperature. On the other hand, 
as expected from the analysis of Fig. [11] the growth of 
domains in BMLJ appears to be more limited. 

Such domains are expected to provide an efficient 
mechanism of slowing down in supercooled liquids.^ To 
address this issue, we evaluate the self intermediate scat- 
tering function for small particles [see Eq. ([3])], according 
to their role in the locally preferred structure at time 
t — tQ. At each time origin to, we partition the small 
particles into c-, n-, and o-particles, as described above. 
The correlation functions F^{k,t), F^{k,t), and F°{k,t) 
are then obtained by performing the average over time 
origins in Eq. ([3]) using only c-, n-, or o-particles respec- 
tively. Relaxation times t^, t„, and Tq are defined as 
in Sec. IIVI The ratio Tc/tq provides a simple measure 
of the slowness of particles inside domains formed by a 
particular locally preferred structure. The temperature 
dependence of Tc/tq is shown in Fig. [T3]for various mix- 
tures at P = 10. The value of Tc/tq tends to 1 at high 
temperature in all systems and increases more markedly 
by decreasing temperature, as the fragility of the sys- 
tem increases. Around T^, we find that relaxation times 



within icosahedral domains differ by roughly an order of 
magnitude from those outside, whereas prismatic struc- 
tures in non-additive mixtures develop a more modest 
separation of time scales. 

The dynamical impact of locally preferred structures is 
assisted, at low temperature, by an increased lifetime of 
such slow domains. To address this point we proceeded 
similarly to Donati et al.^ introducing a single-particle 
function Viit) that equals 1 if particle i is at the center 
of a given Voronoi polyhedron, and if not. Restricting 
our attention to small particles, we computed the auto- 
correlation function— 



N2 



N2 



(8) 



where Tip — the average number of 

small particles at the center of a given polyhedron. We 
estimated the lifetime Tp of a polyhedron from the condi- 
tion <j{Tp) = CT(0)/e. Independent of the polyhedron un- 
der consideration, the normalized autocorrelation func- 
tion a{t)/a{Q) falls rapidly to zero in the normal liquid 
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FIG. 14: Dynamical impact of locally preferred structures, as 
identified by the temperature dependence of the ratio Tc/to 
(main plot) and Tn/ro (insets) at P = 10. See text for def- 
inition of Tc, Tn, and To. Upper plot: AMLJ-A mixtures for 
A = 0.60,0.70,0.73,0.82. Lower plot: BMLJ (filled circles), 
WAHN (open circles) and NiasYer (stars). The dotted line 
drawn at 1 indicates the high temperature limit. 



regime. As the temperature is lowered, polyhedra cor- 
responding to locally preferred structures become more 
long-lived than the others, as expected. Around T^, we 
find that a{t)/a{0) for locally preferred structures de- 
cays to zero within the time scale given by the decay 
of F^{k*,t). In Table [TVl we report the lifetimes Tp 
of some frequent polyhedra found at the lowest equili- 
brated temperatures for WAHN, BMLJ, and NisaYgr . 
In WAHN and BMLJ, the lifetime of polyhedra corre- 
sponding to locally preferred structures is around an or- 
der of magnitude larger than those of other geometries. 
Interestingly, in the case of NiaaYgr , we find that some 
less frequent polyhedra, such as (0,2,8)-polyhedra, have 
a lifetime comparable to that of our putative locally pre- 
ferred structure, suggesting the existence of competing 
structures. We also find that icosahedra tend to have a 
longer lifetime, relative to the typical structural relax- 
ation times, than prismatic structures. 

The relation between fragility and local order, which is 



TABLE IV: Lifetime Tp of most frequent Voronoi polyhedra 
around small particles. Results are obtained from local min- 
ima at the lowest equilibrated temperatures (T ~ Tr). Also 
shown is the ratio Tpjr^, where T2 is the relaxation time ob- 
tained from the condition F^{k'',T) — l/e. 
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apparent from our simulation data, fits rather well into 
the scenario of the frustration-limited domains theory.— 
According to this approach, glass-formation arises from 
the competition of a tendency to form mesoscopic, stable 
domains, characterized by locally preferred structures, 
and a mechanism of frustration, which prevents these 
domains from tiling three dimensional space. Fragility 
turns out to be proportional to the energetic stability of 
such domains and inversely proportional to the strength 
of frustration. At present, the roles of stability and frus- 
tration cannot be clearly disentangled. Nevertheless, the 
following considerations, based on the present work, are 
possible and we hope they could serve as guidelines for 
further theoretical modeling or investigations: (i) In the 
case of additive mixtures, within the explored range of 
size ratio, icosahedral ordering seems to be the most 
prominent structural feature. Results obtained for iso- 
lated Lennard- Jones clustersSS, suggest that the maxi- 
mum of icosahedral ordering around A « 0.84 might be 
related to an enhanced energetic stability of equimolar 
icosahedra, i.e. icosahedra with the same number of large 
and small neighbors. Formation of more stable icosahe- 
dral structures would, in turn, explain the increase of 
fragility with size ratio. A more detailed study of larger 
clusters in the bulk, forming extended regions of icosahe- 
dral coordination,^'^ would probably be required to fur- 
ther clarify this point, (ii) It should be realized that 
frustration in different systems may be of different ori- 
gin. Icosahedral ordering, while being frustrated from 
tiling three-dimensional space for geometrical reasons, al- 
lows the growth of relatively large domains, when com- 
pared to prismatic structures observed in non-additive 
Lennard- Jones mixtures. In non- additive alloys different 
competing locally preferred structures and mismatch in 
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stoichiometry may further increase frustration. 

VI. CONCLUSIONS 

In this work, we have provided a comparative study of 
different supercooled Lennard-Jones mixtures by quench- 
ing the systems at constant pressure. The models an- 
alyzed include the well-studied mixtures BMLJ^i and 
WAHN,"2S. a set of equimolar, additive mixtures with 
varying size ratio, and a model meant the give a real- 
istic structural description of the NiasYer alloy.— These 
systems display a varying degree of fragility, which has 
been rationalized in terms of the properties of some rel- 
evant locally preferred structures. 

Local order has been characterized using a Voronoi 
construction. Employing an effective definition of the 
locally preferred structure of a liquid, as the geometry 
corresponding to the most frequent Voronoi polyhedra, 
we have shown that fragility is related, in the mixtures 
considered in this work, to the rapid growth with tem- 
perature of slow, stable domains characterized by the lo- 
cally preferred structure of the mixture, generalizing pre- 
vious observations on monoatomic bulk systemsrS*^ Such 
a growth with temperature is more rapid, the more fragile 
the mixture. 

Analyzing the set of AMLJ-A mixtures, we found that 
the size ratio A controls the formation of icosahedral or- 
dering in the bulk. Extended regions of icosahedral co- 
ordination are more rapidly formed upon supercooling of 
additive mixtures with moderate size asymmetry. This, 
in turn, leads to a more pronounced super- Arrhenius be- 
havior. These results might also be representative of a 
wider class of mixtures with additive, non-directional in- 
teractions, since icosahedral ordering has been found to 



display a similar trend of variation with size ratio in mod- 
els of bidisperse Cu glasses4^ On the other hand, non- 
additive mixtures, such as BMLJ and NiaaYe? , favor the 
formation of prismatic structures. The growth, upon su- 
percooling, of these structures is milder compared to that 
of icosahedra in additive mixtures. Consequently, non- 
additive mixtures display a less fragile behavior, which 
might be due to the presence of stronger frustration 
mechanisms. From such discussion of our data, it is thus 
tempting to relate fragility to the thermal rate of growth 
of locally preferred structures, whatever their type may 
be. 

In the light of the frustration-limited domains theory;^ 
fragility should originate from the interplay between sta- 
bility of domains formed by locally favored geometries 
and frustration. Our results indicate that further studies 
are needed to assess the relative role of these two factors 
in supercooled mixtures. Deeper investigations on the 
link between fragility and local order, either by study of 
isolated clusters^S or by direct determination of locally 
preferred structures in the bulkiiii^ or by confocal mi- 
croscopy technique)^ will certainly be rewarding. 
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